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We present a comprehensive study of the influence of the geomagnetic field on the energy 
estimation of extensive air showers with a zenith angle smaller than 60°, detected at the 
Pierre Auger Observatory. The geomagnetic field induces an azimuthal modulation of the 
estimated energy of cosmic rays up to the ~ 2% level at large zenith angles. We present a 
method to account for this modulation of the reconstructed energy. We analyse the effect 
of the modulation on large scale anisotropy searches in the arrival direction distributions of 
cosmic rays. At a given energy, the geomagnetic effect is shown to induce a pseudo-dipolar 
pattern at the percent level in the declination distribution that needs to be accounted for. 
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1 Introduction 



High energy cosmic rays generate extensive air showers in the atmosphere. The trajectories of the 
charged particles of the showers are curved in the Earth's magnetic field, resulting in a broadening 
of the spatial distribution of particles in the direction of the Lorentz force. While such effects are 
known to distort the particle densities in a dramatic way at zenith angles larger than ~60° [TJ ® 
[31 0] , they are commonly ignored at smaller zenith angles where the lateral distribution function 
is well described by empirical models of the NKG-type [51 [3] based on a radial symmetry of the 
distribution of particles in the plane perpendicular to the shower axis. 

In this article, we aim to quantify the small changes of the particle densities at ground induced 
by the geomagnetic held for showers with zenith angle smaller than ~60°, focusing on the impacts 
on the energy estimator used at the Pierre Auger Observatory. As long as the magnitude of these 
effects lies well below the statistical uncertainty of the energy reconstruction, it is reasonable to 
neglect them in the framework of the energy spectrum reconstruction. As the strength of the 
geomagnetic held component perpendicular to the arrival direction of the cosmic ray, St, depends 
on both the zenith and the azimuthal angles (#, if) of any incoming shower, these effects are 
expected to break the symmetry of the energy estimator in terms of the azimuthal angle if. Such 
an azimuthal dependence translates into azimuthal modulations of the estimated cosmic ray event 
rate at a given energy. For any observatory located far from the Earth's poles, any genuine large 
scale pattern which depends on the declination translates also into azimuthal modulations of the 
cosmic ray event rate. Thus to perform a large scale anisotropy measurement it is critical to account 
for azimuthal modulations of experimental origin and for those induced by the geomagnetic field, 
as already pointed out in the analysis of the Yakutsk data [7] and the ARGO-YBJ data |S]. Hence, 
this work constitutes an accompanying paper of a search for large scale anisotropics, both in right 
ascension and declination of cosmic rays detected at the Pierre Auger Observatory, the results of 
which will be reported in a forthcoming publication. 

To study the influence of the geomagnetic field on the cosmic ray energy estimator, we make 
use of shower simulations and of the measurements performed with the surface detector array 
of the Pierre Auger Observatory, located in Malargue, Argentina (35.2°S, 69.5°W) at 1400 m 
a.s.l. [9]. The Pierre Auger Observatory is designed to study cosmic rays (CRs) with energies 
above ~ 10 18 eV. The surface detector array consists of 1660 water Cherenkov detectors sensitive 
to the photons and the charged particles of the showers. It is laid out over an area of 3000 km 2 
on a triangular grid and is overlooked by four fluorescence detectors. The energy at which the 
detection efficiency of the surface detector array saturates is ~ 3EeV [TP] . For each event, the 
signals recorded in the stations are fitted to find the signal at 1000 m from the shower core, 5(1000), 
used as a measure of the shower size. The shower size 5(1000) is converted to the value 538 that 
would have been expected had the shower arrived at a zenith angle of 38°. 538 is then converted 
into energy using a calibration curve based on the fluorescence telescope measurements [11) . 

The influence of the geomagnetic field on the spatial distribution of particles for showers with 
zenith angle less than 60° is presented in Section [2J through a toy model aimed at explaining 
the directional dependence of the shower size 5(1000) induced by the geomagnetic field. The 
observation of this effect in the data of the Pierre Auger Observatory is reported in Section [3l 
In Section HI we quantify the size of the 5(1000) distortions with zenith and azimuthal angles by 
means of end-to-end shower simulations, and then present the procedure to convert the shower 
size corrected for the geomagnetic effects into energy using the Constant Intensity Cut method. In 
Section[5j the consequences on large scale anisotropies are discussed, while systematic uncertainties 
associated with the primary mass, the primary energy and the number of muons in showers are 
presented in Section [5] 
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Figure 1: The shower front plane coordinate system [2 [4]: e z is anti-parallel to the shower direction u, 
while e y is parallel to Bt, the projection of the magnetic field B onto the shower plane x-y. (ip, r) are the 
polar coordinates in the shower plane. 

2 Influence of the geomagnetic field on extensive air showers 

The interaction of a primary cosmic ray in the atmosphere produces mostly charged and neutral 
pions, initiating a hadronic cascade. The decay of neutral pions generates the electromagnetic 
component of the shower, while the decay of the charged pions generates the muonic one. Elec- 
trons undergo stronger scattering, so that the electron distribution is only weakly affected by the 
geomagnetic deflections. Muons are produced with a typical energy E^ of a few GeV (increasing 
with the altitude of production). The decay angle between pions and muons is causing only a 
small additional random deflection, as they almost inherit the transverse momentum of their 
parents (a few hundred MeV/c) so that the distance of the muons from the shower core scales as 
the inverse of their energy. While the radial offset of the pions from the shower axis is of the order 
of a few 10 m, it does not contribute significantly to the lateral distribution of the muons observed 
on the ground at distances r > 100 m. Hence, at ground level, the angular spread of the muons 
around the shower axis can be considered as mainly caused by the transverse momentum inherited 
from the parental pions. 

After their production, muons are affected by ionisation and radiative energy losses, decay, 
multiple scattering and geomagnetic deflections. Below 100 GeV, the muon energy loss is mainly 
due to ionisation and is relatively small (amounting to about 2 MeV g _1 cm 2 ), allowing a large 
fraction of muons to reach the ground before decaying. Multiple scattering in the electric field of 
air nuclei randomises the directions of muons to some degree, but the contribution to the total 
angular divergence of the muons from the shower axis remains small up to zenith angles of the 
shower-axis of about 80°. 

Based on these general considerations, we now introduce a simple toy model aimed at un- 
derstanding the main features of the muon density distortions induced by the geomagnetic field. 
We adopt the shower front plane coordinate system depicted in Fig. [1] [2]. In the absence of 
the magnetic field, and neglecting multiple scattering, a relativistic muon of energy ~ cp^ and 
transverse momentum p^ will reach the shower front plane after traveling a distance d at a position 
r from the shower axis given by 

r^d H ?i (1) 
On the other hand, in the presence of the magnetic field, muons suffer additional geomagnetic 
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Figure 2: Magnetic deviations as a function of the distance to the shower axis observed on a simulated 
vertical shower (points). Superimposed are the deviations expected from Eq. © (line). The shaded region 
and the error bars give the corresponding dispersion. 

deflections. We treat the geomagnetic field B in Malargiie as a constant field3, 

B = 24.6^T, L> B = 2.6°, J B = -35.2°, (2) 

Db and Iq being the geomagnetic declination and inclination. The deflection of a relativistic muon 
in the presence of a magnetic field with transverse component Bj> can be approximated with 

ecB T d 2 

^ 2E ' ® 

where e is the elementary electric charge and the sign corresponds to positive/negative charged 
muons. The dependence of the geomagnetic deflections Sx = Sx + = —5x- on the distance to the 
shower axis r — \J x 2 + y 2 is illustrated in Fig. [2] obtained by comparing the position of the same 
muons in the presence or in the absence of the geomagnetic field in a simulated vertical shower of a 
proton at 5 EeV. The deviations expected from the expression for Sx± are also shown in the same 
graph (solid line). It was obtained by inserting muon energy and distance at the production point 
of the simulated muons into Eq. ([3]). It turns out that Eq. © estimates rather well the actual 
deviations, though the distance between the actual and the predicted deviations increases at large 
r. This is mainly because on the one hand d underestimates the actual travel length to a larger 
extent at larger r, while on the other hand the magnetic deviation actually increases while muons 
gradually lose energy during travel. Hence, from the muon density p^(x, y) in the transverse plane 
in the absence of the geomagnetic field, the corresponding density p„(S, y) in the presence of such 
a field can be obtained by making the following Jacobian transformation, in the same way as in 
the framework of very inclined showers [2] , 

Pn(x(x,y),y(x,y)). (4) 

Here, the term "muon density" refers to the time-integrated muon flux through the transverse 
shower front plane associated to the air shower, and the barred coordinates represent the positions 
of the muons in the transverse plane in the presence of the geomagnetic field: 

x = x + Sx±(x,y), 

V = V- (5) 

Mil Malargiie the geomagnetic field has varied by about 1° in direction and 2% in magnitude over 10 years 1121 . 
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Figure 3: Relative changes of Ap M /p M in the transverse shower front plane due to the presence of the 
geomagnetic field, obtained at zenith angle 9 — 60° and azimuthal angle aligned along Db + 180°. 



Since Eq. (Q} induces changes of the shower size 5(1000), it is of particular interest to get an 
approximate relationship between p and p around 1000 m. From Fig. [2 it is apparent that around 
1000 m the mean magnetic deviation is approximately constant over a distance range larger than 
the size of the deviation. When focusing on the changes of density at 1000 m from the shower 
core, it is thus reasonable to neglect the x and y dependence of the deviation 5x±, which allows 
an approximation of the density p M (ir, y) around 1000 m as 

P^(x,V) - Pn+(x- Sx+,y) +p fJ ,_{x- Sx^,y) 

^ P^(x,y) + — — ^r{x,y), (6) 
2 ox 

where we assumed p M _ = p ll+ = p^/2. The two opposite muon charges cancel out the linear 
term in Sx and we see that magnetic effects change the muon density around 1000 m by a factor 
proportional to (5x) 2 oc B\ oc sin 2 (u, b), where u and b = B/|B| denote the unit vectors in the 
shower direction and the magnetic field direction, respectively. This is particularly important with 
regard to the azimuthal behaviour of the effect, as the azimuthal dependence is contained only in 
the B^(9, ip) term. This dependency is therefore a generic expectation outlined by this toy model. 
The model will be verified in Section U by making use of complete simulation of showers. On 
the other hand, the zenith angle dependence relies on other ingredients that we will probe in an 
accurate way in Section 2J such as the altitude distribution of the muon production and the muon 
energy distribution. 

3 Observation of geomagnetic effects in the Pierre Auger 
Observatory data 

To illustrate the differences between p^ and p M described in Eq. (0}, the relative changes Ap^/p^ 
are shown in Fig. |3]in the transverse shower front plane by producing muon maps from simulations 
at zenith angle 9 = 60° and azimuthal angle aligned along Db + 180° in the presence and in the 
absence of the geomagnetic field. A predominant quadrupolar asymmetry at the few percent level 
is visible, corresponding to the separation of positive and negative charges in the direction of the 
Lorentz force. 

This quadrupolar asymmetry is expected to induce to some extent a quadrupolar modulation 
of the surface detector signals as a function of the polar angle on the ground, defined here as the 
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Figure 4: Definition of angle $ with respect to the magnetic East -B mag and the shower core for a given 
shower direction u and a surface detector at r. The azimuthal angle of the magnetic field vector B defines 
the magnetic North AT mag . 
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Figure 5: Average ratio of the true signal in each surface detector with respect to the expected one as 
a function of the polar angle on the ground. Left panel: using simulated showers in the presence (thick 
points) and in the absence (thin points) of the geomagnetic field. Right panel: using real data above 
4 EeV. The solid lines give the fit of a quadrupolar modulation to the corresponding points. 



angle between the axis given by the shower core and the surface detector, and the magnetic East 
y>3 = —D-q = —2.6° (Fig. 0]). The use of this particular angle, instead of the polar angle ip which 
is defined in the shower front plane (see Fig. [T]), allows us to remove dipolar asymmetries in the 
surface detector signals, the origin of which is related to the radial divergence of particles from 
the shower axis. Such asymmetries cancel out in this analysis, due to the isotropic distribution 
of the cosmic rays. To demonstrate the geomagnetic effect, we produced a realistic Monte-Carlo 
simulation using 30 000 isotropically distributed showers (with zenith angles less than 60°) with 
random core positions within the array. The injected primary energies were chosen to be greater 
than 4 EeV (safely excluding angle dependent trigger probability) and distributed according to a 
power law energy spectrum dN/dE oc E" 1 with power index 7 = 2.7, so that this shower library is 
as close as possible to the real data set. To each shower we apply the reconstruction procedure of 
the surface detector, leading to a fit of the lateral distribution function [II] . The lateral distribution 
function parametrizes the signal strength in the shower plane, assuming circular shower symmetry. 
By evaluating the lateral distribution function at the position of the surface detector, we obtain the 
expected signal Scxp- This signal can be compared to the true signal in the surface detector Ssd- 
The ratio between the observed and expected signals as a function of the polar angle on the ground 
in simulated showers is shown in the left panel of Fig. [5j with (thick points) and without (thin 
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points) the geomagnetic field. While a significant quadrupolar modulation with a fixed phase along 
Db and amplitude ~ (1.1 ± 0.2)% is observed when the field is on, no such modulation is observed 
when the field is off (~ (0.1±0.2)%), as expected. In the right panel, the same analysis is performed 
on the real data above 4 EeV, including again about 30 000 showers. A significant modulation of 
~ (1.2 ± 0.2)% is observed, agreeing both in amplitude and phase within the uncertainties with 
the simulations performed in the presence of the geomagnetic field. This provides clear hints of 
the influence of the geomagnetic field in the Auger data. 

Note that this analysis is restricted to surface detectors that are more than 1000 m away from 
the shower core. This cut is motivated by Fig. [3J showing that the quadrupolar amplitude is 
larger at large distances from the shower core. We further require the surface detectors to have 
signals larger than 4 VEMB This cut is a compromise between keeping good statistics and keeping 
trigger effects small. Above 4 VEM the measured amplitude does not depend systematically on 
the signal strength cut. However a cut in the surface detector signals induces a statistical trigger 
bias because showers with upward signal fluctuations will trigger more readily. This explains the 
small discrepancy between real and Monte-Carlo data in terms of the global normalisation in Fig. [5] 
which differs from 1 by ~3%. Cutting at larger signals reduces this discrepancy. 

Most importantly, depending on the incoming direction, the quadrupolar asymmetry is also 
expected to affect the shower size 5(1000) and thus the energy estimator as qualitatively described 
in Eq. ©. Consequently, these effects are expected to modulate the estimated cosmic ray event 
rate at a given energy as a function of the incoming direction, and in particular to generate a 
North/South asymmetry in the azimuthal distribution^. Such an asymmetry is also expected in 
the case of a genuine large scale modulation of the flux of cosmic rays. However related analyses 
of the azimuthal distribution are out of the scope of this paper, and we restrict ourselves in the 
rest of this article to present a comprehensive study of the geomagnetic distortions of the energy 
estimator. This will allow us to apply the corresponding corrections in a forthcoming publication 
aimed at searching for large scale anisotropies. 

4 Geomagnetic distortions of the energy estimator 
4.1 Geomagnetic distortions of the shower size S'(1000) 

The toy model presented in Section [2] allows us to understand the main features of the influence 
of the geomagnetic field on the muonic component of extensive air showers. To get an accurate 
estimation of the distortions induced by the field on the shower size S'(IOOO) as a function of both 
the zenith and the azimuthal angles, we present here the results obtained by means of end-to- 
end simulations of proton-initiated showers generated with the AIRES program [14] and with the 
hadronic interaction model QGSJET [TS]. We have checked that the results obtained with the 
CORSIKA program [TB] are compatible. We consider a fixed energy E = 5 EeV and seven fixed 
zenith angles between 8 = 0° and 8 = 60°. The dependency of the effect in terms of the primary 
mass and of the number of muons in showers as well as its evolution with energy are sources 
of systematic uncertainties. The influence of such systematics will be quantified in Section [6l 
Within our convention for the azimuthal angle, the azimuthal direction of the magnetic North is 

= 90° - D B = 87.4°. The zenith direction of the field is 8 B = 90° - |J B | = 54.8°. 

To verify the predicted behaviour of the shower size shift in terms of B^, we first show the 
results of the simulations of 1000 showers at a zenith angle 8 = 8b and for two distinct azimuthal 
angles cp = ip^ and ip = (p^ +90°. Each shower is then thrown 10 times at the surface detector array 
with random core positions and reconstructed using exactly the same reconstruction procedure as 

2 VEM - Vertical Equivalent Muon - is the average charge corresponding to the Cherenkov light produced by a 
vertical and central through-going muon in the surface detector. It is the unit used in the evaluation of the signal 
recorded by the detectors 

3 The convention we use for the azimuthal angle tp is to define it relative to the East direction, counterclockwise. 
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Figure 6: Distributions of shower size ,5(1000) obtained by simulating showers at zenith angle 9 — 9b and 
azimuthal angle ip^ (left) and ifB +90° (right). Thick histogram: no magnetic field. Dotted histogram: 
real magnetic field in Malargiie. Dashed histogram: twice the real magnetic field in Malargiie. 
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Figure 7: A5 , (1000)/5'(1000) (in %) as a function of the azimuthal angle tp, at zenith angle 9 = 9b for 
two different field strengths. Points are obtained by Monte Carlo shower simulation, lines are the expected 
behavior (see Section [2j|. 



the one applied to real data. For this specific zenith angle 6*b, no shift is expected in the North 
direction ip^ as the transverse component of the magnetic field is zero. This is indeed the case as 
illustrated in the left panel of Fig. [BJ showing the distribution of reconstructed S'(IOOO) for three 
different configurations of the magnetic field: no field, real field in Malargiie, and twice the real field 
in Malargiie. It can be seen that on average all histograms are - within the statistical uncertainties 
on the average - centered on the same value. In the right panel of Fig. |6] we repeat the same 
analysis with the showers generated in the direction ip^ + 90° . Since the transverse component of 
the field is now different from zero, a clear relative shift in terms of A5(1000)/5(1000) is observed 
between the three distributions: the shift is ~ 1.6% between the configurations with and without 
the field, leading to a discrimination with a significance of ~ 5.5 a, while the shift is ~ 6% between 
the configurations with twice the real field and without the field leading to a discrimination with a 
significance of ~ 20 a. It can be noticed that the strength of the shift is thus in overall agreement 
with the expected scaling B\. 

For the zenith angle 9 — 6b, in Fig. [7] we show the shift of the mean 5(1000) obtained by 
simulating 1000 showers in the same way as previously for eight different values of the azimuth 
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Figure 8: G(9) = A5'(1000)/S(1000)/ sin 2 (u, b) as a function of the zenith angle 9. 



angle. Again, the results are displayed for configurations with the real field (bottom) and with twice 
the real field (top). The expected behaviours in terms of A5(1000)/5(1000) = G{9 B ) sin 2 (ujb) 
are shown by the continuous curves, where the normalisation factor G is tuned by hand. Clearly, 
the shape of the curves agrees remarkably well with the Monte Carlo data within the uncertainties. 
Hence, this study supports the claim that the azimuthal dependence of the shift in 5(1000) induced 
by the magnetic field is proportional to B^(9, ip), in agreement with the expectations provided by 
general considerations expressed in the previous section on the muonic component of the showers. 

The B\ term encompassing the overall azimuthal dependence at each zenith angle, the remain- 
ing shift G(9) = A5(1000)/5(1000)/sin 2 (u7b) depends on the zenith angle through the altitude 
distribution of the muon production, the muon energy distribution, and the weight of the muonic 
contribution to the shower size 5(1000). Repeating the simulations at different zenith angles, we 
plot G as a function of the zenith angle in Fig. [5] Due to the increased travel lengths of the muons 
and due to their larger relative contribution to 5(1000) at high zenith angles, the value of G rises 
rapidly for angles above ~ 40°. The superimposed curve is an empirical fit, allowing us to get the 
following parametrisation of the shower size distortions induced by the geomagnetic field, 

A5(1000), a 1n _ 3 _-2.8 



5(1000) 



(0, ip) = 4.2 ■ 10-* cos 9 sin^(u, b). (7) 



4.2 Prom shower size to energy 

At the Pierre Auger Observatory, the shower size 5(1000) is converted into energy E using a two- 
step procedure [TT|. First, the evolution of 5(1000) with zenith angle arising from the attenuation 
of the shower with increasing atmospheric thickness is quantified by applying the Constant Inten- 
sity Cut (CIC) method that is based on the (at least approximate) isotropy of incoming cosmic 
rays. The CIC relates relates 5(1000) in vertical and inclined showers through a line of equal 
intensity in spectra at different zenith angles. This allows us to correct the value of 5(1000) for 
attenuation by computing its value had the shower arrived from a fixed zenith angle, here 38 de- 
grees (corresponding to the median of the angular distribution of events for energies greater than 
3 EeV). This zenith angle independent estimator 538 is defined as 538 = S (1000) /CIC (9). The 
calibration of 53§ with energy E is then achieved using a relation of the form E = AS^ 8 , where 
A = 1.49 ± 0.06(stat)±0.12(syst) and B = 1.08 ± 0.01(stat)±0.04(syst) were estimated from the 
correlation between 538 and E in a subset of high quality "hybrid events" measured simultaneously 
by the surface detector (SD) and the fluorescence detector (FD) [TT]. In such a sample, 538 and E 
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Figure 9: Relative differences AN/N as a function of the declination, for 2 different values of # ma x- 

are independently measured, with S^s from the SD and E from the FD. 

This two-step procedure has an important consequence on the implementation of the energy 
corrections for the geomagnetic effects. The CIC curve is constructed assuming that the shower 
size estimator S'(IOOO) does not depend on the azimuthal angle. The induced azimuthal variation 
of S'(IOOO) due to the geomagnetic effect is thus averaged while the zenith angle dependence of the 
geomagnetic effects is absorbed when the CIC is implemented. To illustrate this in a simplified 
way, consider the case in which the magnetic field were directed along the zenith direction {i.e. in 
the case of a virtual Observatory located at the Southern magnetic pole) so that the transverse 
component of the magnetic field would not depend on the azimuthal direction of any incoming 
shower. Then the shift in ,5(1000) would depend only on the zenith angle in such a way that the 
Constant Intensity Cut method would by construction absorb the shift induced by G(9) into the 
empirical CIC (9) curve, while the empirical relationship E = AS^ g would calibrate S^s into energy 
with no need for any additional corrections. 

This leads us to implement the energy corrections for geomagnetic effects, relating the energy 
Eq reconstructed ignoring the geomagnetic effects to the corrected energy E by 



with 



where 



E 



A(8,<p) = G(6) 



E 



(1 + A(0,<p)y 



sin (u,b) — (sin (u, b) 



(8) 



(9) 



denotes the average over ip and where B is one of the parameters used in the S^s to 
E conversion described above. This expression implies that energies are under- estimated preferen- 
tially for showers coming from the northern directions of the array, while they are over-estimated 
for showers coming from the southern directions, the size of the effect increasing with the zenith 
angle. 



5 Consequences for large scale anisotropy searches 
5.1 Impact on the estimated event rate 

To provide an illustration of the impact of the energy corrections for geomagnetic effects, we 
calculate here, as a function of declination 5, the deviation of the event rate Nq(6), measured if we 
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were not to implement the corrections of the energy estimator by Eq. ©, to the event rate N(6) 
expected from an isotropic background distribution. 

The "canonical exposure" [IT] holds for a full-time operation of the surface detector array above 
the energy at which the detection efficiency is saturated over the considered zenith range. In such 
a case, the directional detection efficiency is simply proportional to cos 6?, 

u(6) oc cos(0) H(6 - ma *) (10) 

where H is the Heaviside function and # max is the maximal zenith angle considered. The zenith 
angle is related to the declination 5 and the right ascension a through 

cos 9 = sin ^ s ite sin S + cos 4ite cos <5 cos a (11) 

where 4itc is the Earth's latitude of the Observatory. The event rate at a given declination 6 and 
above an energy threshold E t h is obtained by integrating in energy and right ascension a, 

NW«r*E r^(0) dNie ^ E) (12) 



dE 



Note that at lower energies this integral acquires an additional energy and angle dependent detec- 
tion efficiency term t{E, 9, 4>). Hereafter we assume that the cosmic ray spectrum is a power law, 
i.e. dN/dE cx E~ 7 . From Eq. ((5J it follows that if the effect of the geomagnetic field were not 
accounted for, the measured energy spectrum would have a directional modulation given by 

^tx[l + A(0, V )] B ^Ep. (13) 

This leads to the following measured event rate above a given uncorrected energy £? t h» 

N {S)oc dE daH(cos9-cos9 max )cos9[l + A(9, l f)} B{ ~ f ~ 1) E~~ f , (14) 

J Eth J 

where ip is related to a and 6 through 

sin<Jcos4ite - cos£cosasin4itc 

tanw= — . (15) 

cos o sin a 



The event rate N (6) as a function of declination is then calculated using Eq. flT3l in Eq. (fl! 
The relative difference AN/N is shown in Fig. |H] as a function of the declination, with spectral 
index 7 = 2.7. The energy over-estimation (under-estimation) of events coming preferentially from 
the Southern (Northern) azimuthal directions, as described in Eq. JSJ, leads to an effective excess 
(deficit) of the event rate for S < —20° (S > —20°), with an amplitude of ~ 2% when considering 
= 60°. It is worth noting that this amplitude is reduced to within 1% when considering 
50°, as shown by the dotted line. 
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5.2 Impact on dipolar modulation searches 

The pattern displayed in Fig. [9] roughly imitates a dipole with an amplitude at the percent level. 
To evaluate precisely the impact of this pattern on the assessment of a dipole moment in the 
reconstructed arrival directions and to probe the statistics needed for the sensitivity to such a 
spurious pattern, we apply the multipolar reconstruction adapted to the case of a partial sky 
coverage p"8] to mock data sets by limiting the maximum bound of the expansion L max to 1 (pure 
dipolar reconstruction). Since the distortions are axisymmetric around the axis defined by the 
North and South celestial poles, only the multipolar coefficient related to this particular axis is 
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Figure 10: Dipolar reconstruction of arrival directions of mock data sets with event rates distorted by the 
geomagnetic effects. Left: distributions of amplitudes. Right: distributions of declinations. The smooth 
lines give the expected distribution in the case of isotropy. 



expected to be affected (here: ain). Consequently, this particular coefficient has impacts on both 
the amplitude of the reconstructed dipole and its direction with respect to the axis defined by the 
North and South celestial poles (the technical details of relating the estimation of the multipolar 
coefficients to the spherical coordinates of a dipole are given in the Appendix). 

To simulate the directional distortions induced by Eq. ([SJ, each mock data set is drawn from 
the event rate Nq(S) corresponding to the uncorrected energies, and is reconstructed using the 
canonical exposure in Eq. (|10j) . The results of this procedure applied to 1000 samples are shown 
in Fig. [TU] In the left panel, the distribution of the reconstructed amplitudes r using N = 300 000 
events is shown by the dotted histogram. It clearly deviates from the expected isotropic distribution 
displayed as the dotted curve which corresponds to (see Appendix) 



Pn(r) = erfi(^HZ^=) exp ( - £) , (16) 

where erfi(z) = erf (iz)/i, and where the width parameters a and a z can be calculated from 
the exposure function [T^- With the particular exposure function used here, it turns out that 
a ~ 1.02-y/3/iV and a z ~ 1.59^/3/ N. This allows us to estimate the spurious dipolar amplitudc0 
to be of the order of the mean of the dotted histogram, about ~ 1.9%. Consequently, we can 
estimate that the spurious effect becomes predominant as soon as the mean noise amplitude (r) 
deduced from Eq. (fTB"j) is of the order of 1.9%, 



(r) = K+ ° 2 ™^ l _-/l^ , 1.9%. (17) 

This translates into the condition N ~ 32 000 (solid histogram). Using such a number of events, 
the bias induced on the amplitude reconstruction is illustrated in the same graph by the longer 
tail of the full histogram with respect to the expected one, and is even more evident in the right 
panel of Fig. [TU1 showing the distribution of the reconstructed declination direction of the dipole 
which already deviates to a large extent from the expected distribution. 



4 Due to the partial sky exposure considered here, the estimate of the dipolar amplitude is biased by the higher 
multipolar orders needed to fully describe AN/N shown in Fig. 1101 18]. The aim of this calculation is only to 
provide a quantitative illustration of the spurious measurement which would be performed due to the geomagnetic 
effects when reconstructing a pure dipolar pattern. 
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Figure 11: Relative differences AN/N as a function of the declination, for different primary masses, 
different primary energies, different hadronic models and for increased number of muons in showers. 



6 Systematic uncertainties 

The parametrisation of G(9) in Eq. ([7]) was obtained by means of simulations of proton showers 
at a fixed energy. The height of the first interaction influences the production altitude of muons 
detected at 1000 m from the shower core at the ground level. Moreover, as muons are produced at 
the end of the hadronic cascade, when the energy of the charged mesons is diminished so much that 
their decay length becomes smaller than their interaction length (which is inversely proportional 
to the air density), the energy distribution of muons is also affected by the height of the first 
interaction. Because the air density is lower in the upper atmosphere, this mechanism results in 
an increase of the energy of muons. The muonic contribution to S'(IOOO) depends also on both the 
primary mass and primary energy. For all these reasons, the parametrisation of G{9) is expected 
to depend on both the primary mass and primary energy. 

To probe these influences, we repeat the same chain of end-to-end simulations using proton 
showers at energies of 50 EeV and iron showers at 5 EeV. Results in terms of the distortions of 
the observed event rate N(8) are shown in Fig. [TTJ We also display in the same graph the results 
obtained using the hadronic interaction model QGSJETII [19] , The differences with respect to the 
reference model are small, so that the consequences on large scale anisotropy searches presented 
in Section [5] remain unchanged within the statistics available at the Pierre Auger Observatory. 

In addition, there are discrepancies in the hadronic interaction model predictions regarding 
the number of muons in shower simulations and what is found in our data |20j . Higher number of 
muons influences the weight of the muonic contribution to S'(IOOO). The consequences of increasing 
the number of muons by a factor of 2 on the distortions of the observed event rate are also shown 
in Fig. [TTJ As the muonic contribution to 5(1000) is already large at high zenith angles in the 
reference model, this increase of the number of muons does not lead to large differences. 



7 Conclusion 

In this work, we have identified and quantified a systematic uncertainty affecting the energy deter- 
mination of cosmic rays detected by the surface detector array of the Pierre Auger Observatory. 
This systematic uncertainty, induced by the influence of the geomagnetic field on the shower devel- 
opment, has a strength which depends on both the zenith and the azimuthal angles. Consequently, 
we have shown that it induces distortions of the estimated cosmic ray event rate at a given energy 
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at the percent level in both the azimuthal and the declination distributions, the latter of which 
mimics an almost dipolar pattern. 

We have also shown that the induced distortions are already at the level of the statistical 
uncertainties for a number of events N ~ 32 000 (we note that the full Auger surface detector 
array collects about 6500 events per year with energies above 3 EeV). Accounting for these effects 
is thus essential with regard to the correct interpretation of large scale anisotropy measurements 
taking explicitly profit from the declination distribution. 
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Appendix 

The p.d.f. of the first harmonic amplitude for a data set of N points drawn at random over a 
circle is known to be the Rayleigh distribution. In this appendix, we generalise this distribution 
to the case of N points being drawn at random on the sphere over the exposure lj(6) of the 
Pierre Auger Observatory. Assuming the underlying arrival direction distribution to be of the 
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form $(a, S) = $o(l + D • u), the components of the dipolar vector D are related to the multipolar 
coefficients through 



Oil 

a oo 



ai-i 



fflio 



(18) 



Denoting by x, y, z the estimates of D x , D y , D z , the joint p.d.f. px,Y,z(x, y, z) can be factorised in 
the limit of large number of events in terms of three centered Gaussian distributions N(0,a), 



Px,Y,z(x,y,z) = p x (x)p Y {y)pz(z) = N(0,a x )N(0,a y )N(0,a z ), 



(19) 



where the standard deviation parameters can be calculated from the exposure function [18] . With 
the particular exposure function used here, it turns out that numerical integrations lead to a ~ 
1.02 and a z ~ 1.59^/3/iV. The joint p.d.f. p_r,a. j 4(^, <5, a) expressing the dipole components 
in spherical coordinates is obtained from Eq. (|19[) by performing the Jacobian transformation 



PR,A,A{r,S,a) = 



d(x,y,z) 



d{r, S, a 
r 2 cos 8 
(2tt) 3 /W ; 



Px,y,z(a;(r, 5, a), y(r, 6, a), z(r, 6, a)) 



■ exp 



r 2 cos 2 6 r 2 sin 2 5 



2a 2 



2a 2 



(20) 



From this joint p.d.f., the p.d.f. of the dipole amplitude (declination) is finally obtained by 
marginalising over the other variables, yielding 



Pr{t) 



Pa(S) = 



^2 



: erfi 



V2 



exp 



r 

2^ 



cos S 



2 (a 2 cos 2 5 + a 2 sin 2 <5) 3 / 2 



(21) 



Finally, one can derive from pn quantities of interest, such as the expected mean noise (r), the 
RMS oy and the probability of obtaining an amplitude greater than r: 



- I rr, 

7T 



a 2 arctanh( y/l — cr 2 /a 2 ] 



2a 2 



<jl-(r)\ 



f r \ a ( J a 2 - o 2 \ ( r 2 \ 

Prob(> r) = erfc — — — H erfi ^ r exp — — - 

\V2aJ V^F 3 ^ V V2aa z ) P \ 2a 2 ) 



(22) 

(23) 
(24) 



which are the equivalent to the well known Rayleigh formulas (r) = Ww/N, a r = v/(4 — ir)/N and 
Prob(> r) = cxp(— Nr 2 /A) when dealing with N points drawn at random over a circle [21] . 
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